# R daily, monthly returns: plot, distribution
# VR test for white noise
# resampling distribution for VR test

library("tseries")
library("zoo")
library("moments")

####################################
###  script for basic statistics ###
####################################
basic_stats <- function(ts0) {
  ts <- as.numeric(ts0)
  hh <- c("Basic statistics:    "," Nobs =", formatC(length(ts), width=9) )
  print.noquote(hh)
  hh <- c("Mean=    ", formatC( mean(ts), digits = 5, width= 8, format = "fg") , 
          "StDev=   ", formatC(  sd(ts) , digits = 5, width= 8, format = "fg")  )
  print.noquote(hh )                  
  hh <- c("Skewness=",formatC( skewness(ts), digits = 5, width= 8, format = "fg") ,
          "Kurtosis=",formatC( kurtosis(ts), digits = 5, width= 8, format = "fg")  )
  print.noquote(hh) 
  jarque.bera.test(ts)
}  
#####################################

start0 <- "1998-10-01"
end0   <- "2004-10-01"  

#######################################
### compare two indices, daily data ###
#######################################
# Read data
require(quantmod) # Finanzdaten download

dow <- new.env()
dax <- new.env()

getSymbols("^DJA", env = dow, from = as.Date("1998-10-01"), to = as.Date("2004-10-01"))
getSymbols("^GDAXI", env = dax, from = as.Date("1998-10-01"), to = as.Date("2004-10-01"))
# read in series are of xts type,  plot is different for xts and ts

x.d <-as.zoo( dow$DJA[,4]);       y.d <- as.zoo(dax$GDAXI[,4]); 

## Plot series, returns 
plot(merge(x.d,y.d), main = "Index: ^DJA, ^GDAXI, daily")
#plot(x.d, y.d) 
r.xd  <- diff(log(x.d));         # return US stock index
r.yd  <- diff(log(y.d));         # return other index
plot( merge(r.xd, r.yd))         # when US is coughing, EU has fever

## Distribution of returns, 
plot(density(r.xd), main="r-^DJA ... black, r-^GDAXI ... blue")
lines(density(r.yd), col="blue")  # 2nd plot printed in 1st 

# basic statistics
basic_stats(r.xd)
basic_stats(r.yd)

#########################################
### compare two indices, monthly data ###
#########################################
### monthly data are the end-of-month value:  index is in days!
x.m <- x.d[!duplicated(as.yearmon(time(x.d)), fromLast = TRUE)]
y.m <- y.d[!duplicated(as.yearmon(time(y.d)), fromLast = TRUE)]

#plot(x.m);  plot(y.m);                  
# merge works on daily bases 
# end-of-month-day different in both indices
xm <- ts(as.numeric(x.m), start = c(1998, 10), frequency = 12 )
ym <- ts(as.numeric(y.m), start = c(1998, 10), frequency = 12 )
# merge does not work for ts
x.mz <- as.zoo(xm)
y.mz <- as.zoo(ym)
plot(merge(x.mz,y.mz), main = "Index: ^DJA, ^GDAXI, monthly")
#plot(x.mz, y.mz) 

r.xmz  <- diff(log(x.mz));         # return US stock index
r.ymz  <- diff(log(y.mz));         # return other index

#plot(r.xmz); plot(r.ymz);
plot( merge(r.xmz, r.ymz), main = "Return: ^DJA, ^GDAXI, monthly")      

## Distribution of returns, 
plot(density(r.xmz), main="r-^DJA ... black, r-^GDAXI ... blue")
lines(density(r.ymz), col="blue")        # 2nd plot printed in 1st 

# basic statistics
basic_stats(r.xmz)
basic_stats(r.ymz)

###########################
### Variance Ratio test ###    # based only on daily data 
###########################    # a month is assumed to have 22 trading days
p <- as.vector(log(x.d));       

m <- 22;
M <- as.integer( (length(p)-1) / m ) 
TT <- m*M
p <- p[1:(TT+1)];

## calucation of the sample statistic
r <- diff(p);
#r <- rnorm(TT)           # for checking whether the asy distr is correct
mu <- mean(r)

rm <- c(1:M)*0
for (i in 1:M) { j2 <- m*i; j1 <- m*(i-1)+1;  sumr <- 0;
                 for (ii in j1:j2) {  sumr <- sumr + r[ii]; }
                 rm[i] <- sumr }

s2a <- var(r)*(TT-1)/TT               # var() ... variance corrected for df's 
mcrm <- rm - m*mu                     # mean corrected rm
s2b <- ( var(mcrm)*(M-1)/M )/m
VRm <- s2b/s2a -1           # asy distributed as N(0,2(m-1)/TT)
VRm                        

## resampling of r  ( mu = mean(r), s2a are invariant )
## small sample distribution of VR is simulated
nsim <- 1000
o <- c(1:TT)
VRms <- c(1:nsim)*0

set.seed(1234)
for (is in 1:nsim) 
{
  o   <- sample(o) 
  rs  <- r[o]  
  rms <- c(1:M)*0
  for (i in 1:M) { j2 <- m*i; j1 <- m*(i-1)+1; sumr <- 0; 
                   for (ii in j1:j2) {  sumr <- sumr + rs[ii]; }
                   rms[i] <- sumr }
  
  mcrms <- rms - m*mu                     # mean corrected rm
  s2bs <- ( var(mcrms)*(M-1)/M )/m
  VRms[is] <- s2bs/s2a -1 
}

VRm
quantile(VRms, c(0.05, 0.95) ) 

## comparison, asy and resampling distribution
plot(function(x) dnorm( x,0,(2*(m-1)/TT)^(0.5) ), -1, 1, main="empirical density ... blue ") 
lines(density(VRms), col="blue") 